
library(MASS)
library(Zelig)
library(ZeligChoice)
library(ZeligMultilevel)
library(survival)
library(mlogit)

data<-read.csv("FullData.csv")
nobreak<-read.csv("NoBreak.csv")
water<-read.csv("Water.csv")
survcoop<-read.csv("SurvivalCoop.csv")
survnon<-read.csv("SurvivalNon.csv")
inflow<-read.csv("Inflow.csv")

##########################
##						##
## 		Water plots		##
##						##
##########################

inflow<-subset(inflow,is.na(inflow$Date)==0) #drop extra observations
inflow$Toktogul<-as.numeric(as.character(inflow$Toktogul))
inflow$Nurek<-as.numeric(as.character(inflow$Nurek))


#monthly subsets
jan<-subset(inflow,inflow$Jan==1)
feb<-subset(inflow,inflow$Feb==1)
mar<-subset(inflow,inflow$Mar==1)
apr<-subset(inflow,inflow$Apr==1)
may<-subset(inflow,inflow$May==1)
jun<-subset(inflow,inflow$Jun==1)
jul<-subset(inflow,inflow$Jul==1)
aug<-subset(inflow,inflow$Aug==1)
sep<-subset(inflow,inflow$Sep==1)
oct<-subset(inflow,inflow$Oct==1)
nov<-subset(inflow,inflow$Nov==1)
dec<-subset(inflow,inflow$Dec==1)

#generate monthly averages and standard deviations

syrdarya<-c(mean(jan$Toktogul),mean(feb$Toktogul),mean(mar$Toktogul),mean(apr$Toktogul),mean(may$Toktogul),mean(jun$Toktogul),mean(jul$Toktogul),mean(aug$Toktogul),mean(sep$Toktogul),mean(oct$Toktogul),mean(nov$Toktogul),mean(dec$Toktogul))

syrdarya_sd<-c(sd(jan$Toktogul),sd(feb$Toktogul),sd(mar$Toktogul),sd(apr$Toktogul),sd(may$Toktogul),sd(jun$Toktogul),sd(jul$Toktogul),sd(aug$Toktogul),sd(sep$Toktogul),sd(oct$Toktogul),sd(nov$Toktogul),sd(dec$Toktogul))

syrdarya2<-syrdarya+syrdarya_sd
syrdarya1<-syrdarya-syrdarya_sd

amudarya<-c(mean(jan$Nurek),mean(feb$Nurek),mean(mar$Nurek),mean(apr$Nurek),mean(may$Nurek),mean(jun$Nurek),mean(jul$Nurek),mean(aug$Nurek),mean(sep$Nurek),mean(oct$Nurek),mean(nov$Nurek),mean(dec$Nurek))

amudarya_sd<-c(sd(jan$Nurek),sd(feb$Nurek),sd(mar$Nurek),sd(apr$Nurek),sd(may$Nurek),sd(jun$Nurek),sd(jul$Nurek),sd(aug$Nurek),sd(sep$Nurek),sd(oct$Nurek),sd(nov$Nurek),sd(dec$Nurek))

amudarya2<-amudarya+amudarya_sd
amudarya1<-amudarya-amudarya_sd

#create the deviation variables
inflow$mdevsd_s<-ifelse(inflow$Jan==1,(inflow$Toktogul-mean(jan$Toktogul))/sd(jan$Toktogul),0)
inflow$mdevsd_s<-ifelse(inflow$Feb==1,(inflow$Toktogul-mean(feb$Toktogul))/sd(feb$Toktogul),inflow$mdevsd_s)
inflow$mdevsd_s<-ifelse(inflow$Mar==1,(inflow$Toktogul-mean(mar$Toktogul))/sd(mar$Toktogul),inflow$mdevsd_s)
inflow$mdevsd_s<-ifelse(inflow$Apr==1,(inflow$Toktogul-mean(apr$Toktogul))/sd(apr$Toktogul),inflow$mdevsd_s)
inflow$mdevsd_s<-ifelse(inflow$May==1,(inflow$Toktogul-mean(may$Toktogul))/sd(may$Toktogul),inflow$mdevsd_s)
inflow$mdevsd_s<-ifelse(inflow$Jun==1,(inflow$Toktogul-mean(jun$Toktogul))/sd(jun$Toktogul),inflow$mdevsd_s)
inflow$mdevsd_s<-ifelse(inflow$Jul==1,(inflow$Toktogul-mean(jul$Toktogul))/sd(jul$Toktogul),inflow$mdevsd_s)
inflow$mdevsd_s<-ifelse(inflow$Aug==1,(inflow$Toktogul-mean(aug$Toktogul))/sd(aug$Toktogul),inflow$mdevsd_s)
inflow$mdevsd_s<-ifelse(inflow$Sep==1,(inflow$Toktogul-mean(sep$Toktogul))/sd(sep$Toktogul),inflow$mdevsd_s)
inflow$mdevsd_s<-ifelse(inflow$Oct==1,(inflow$Toktogul-mean(oct$Toktogul))/sd(oct$Toktogul),inflow$mdevsd_s)
inflow$mdevsd_s<-ifelse(inflow$Nov==1,(inflow$Toktogul-mean(nov$Toktogul))/sd(nov$Toktogul),inflow$mdevsd_s)
inflow$mdevsd_s<-ifelse(inflow$Dec==1,(inflow$Toktogul-mean(dec$Toktogul))/sd(dec$Toktogul),inflow$mdevsd_s)

inflow$mdevsd_a<-ifelse(inflow$Jan==1,(inflow$Nurek-mean(jan$Nurek))/sd(jan$Nurek),0)
inflow$mdevsd_a<-ifelse(inflow$Feb==1,(inflow$Nurek-mean(feb$Nurek))/sd(feb$Nurek),inflow$mdevsd_a)
inflow$mdevsd_a<-ifelse(inflow$Mar==1,(inflow$Nurek-mean(mar$Nurek))/sd(mar$Nurek),inflow$mdevsd_a)
inflow$mdevsd_a<-ifelse(inflow$Apr==1,(inflow$Nurek-mean(apr$Nurek))/sd(apr$Nurek),inflow$mdevsd_a)
inflow$mdevsd_a<-ifelse(inflow$May==1,(inflow$Nurek-mean(may$Nurek))/sd(may$Nurek),inflow$mdevsd_a)
inflow$mdevsd_a<-ifelse(inflow$Jun==1,(inflow$Nurek-mean(jun$Nurek))/sd(jun$Nurek),inflow$mdevsd_a)
inflow$mdevsd_a<-ifelse(inflow$Jul==1,(inflow$Nurek-mean(jul$Nurek))/sd(jul$Nurek),inflow$mdevsd_a)
inflow$mdevsd_a<-ifelse(inflow$Aug==1,(inflow$Nurek-mean(aug$Nurek))/sd(aug$Nurek),inflow$mdevsd_a)
inflow$mdevsd_a<-ifelse(inflow$Sep==1,(inflow$Nurek-mean(sep$Nurek))/sd(sep$Nurek),inflow$mdevsd_a)
inflow$mdevsd_a<-ifelse(inflow$Oct==1,(inflow$Nurek-mean(oct$Nurek))/sd(oct$Nurek),inflow$mdevsd_a)
inflow$mdevsd_a<-ifelse(inflow$Nov==1,(inflow$Nurek-mean(nov$Nurek))/sd(nov$Nurek),inflow$mdevsd_a)
inflow$mdevsd_a<-ifelse(inflow$Dec==1,(inflow$Nurek-mean(dec$Nurek))/sd(dec$Nurek),inflow$mdevsd_a)

inflow1<-subset(inflow,Date>=2000)
inflow1$date<-seq(2000,2011-(1/12),1/12)

pdf(file="water.pdf")
par(mfrow=c(3,1))

plot(syrdarya,type="l",ylim=c(0,2000), sub="Syr Darya",ylab="flow",xlab="month",main="Water Flow Monthly Averages, 1992-2010")
lines(syrdarya1,lty="dashed")
lines(syrdarya2,lty="dashed")
legend("topright",c("mean flow","one standard deviation"),lty=c(1,2))
plot(amudarya,type="l",ylim=c(0,2000),sub="Amu Darya",ylab="flow",xlab="month")
lines(amudarya1,lty="dashed")
lines(amudarya2,lty="dashed")
plot(water1$date,water1$mdevsd_s,ylim=c(-3,3),type="l",col="black",ylab="deviation",xlab="year",main="Deviation from Monthly Averages, 2000-2010")
lines(water1$date,water1$mdevsd_a,ylim=c(-3,3),col="black",lty="longdash")
abline(h=0)
legend("bottomright",c("Syr","Amu"),col=c("black","black"),lty=c(1,5))
dev.off()


##############################
## 		  Full-basin: 		##
##		logit analyses  	##
##			 and			##
##			plots			##
##############################


coop1<-zelig(Coop~mdevsd, data=data, model="logit.gee",id="dyad")
coop2<-zelig(Coop~mdevsd+hist_coop12 +grow+secondhalf+exchange + syr, data=data,model="logit.gee",id="dyad")
coop3<-clogit(Coop~mdevsd + strata(dyad), data=data)
coop4<-clogit(Coop~mdevsd+hist_coop12+grow+secondhalf+strata(dyad),data=data)

non1<-zelig(Non~mdevsd, data=data, model="logit.gee",id="dyad")
non2<-zelig(Non~mdevsd +hist_non12+grow+secondhalf +exchange + syr, data=data,model="logit.gee",id="dyad")
non3<-clogit(Non~mdevsd + strata(dyad), data=data)
non4<-clogit(Non~mdevsd+hist_non12+grow+secondhalf + strata(dyad),data=data)


ws<-seq(min(data$mdevsd,na.rm=TRUE),max(data$mdevsd,na.rm=TRUE),0.02)
gr<-mean(data$grow,na.rm=TRUE)
sh<-mean(data$secondhalf,na.rm=TRUE)
ex<-mean(data$exchange,na.rm=TRUE)
sy<-mean(data$syr,na.rm=TRUE)
en<-mean(data$syr,na.rm=TRUE)
hiscoop<-mean(data$hist_coop12,na.rm=TRUE)
hisnon<-mean(data$hist_non12,na.rm=TRUE)

z.out<-coop2
x.out<-setx(z.out,mdevsd=ws,hist_coop12=hiscoop, grow=gr,secondhalf=sh,exchange=ex,syr=sy,energy_diff=en)
s.out<-sim(z.out,x=x.out,num=50000)
pred1<-s.out[[1]][["stats"]][[1]]
k<-length(ws)
	for(i in c(2:k)){
		x<-s.out[[i]][["stats"]][[1]]
		pred1<-rbind(pred1,x)
	}
pred1<-data.frame(pred1,ws)

z.out<-non2
x.out<-setx(z.out,mdevsd=ws, hist_non12=hisnon, grow=gr,secondhalf=sh,exchange=ex,syr=sy,energy_diff=en)
s.out<-sim(z.out,x=x.out,num=50000)
pred2<-s.out[[1]][["stats"]][[1]]
k<-length(ws)
	for(i in c(2:k)){
		x<-s.out[[i]][["stats"]][[1]]
		pred2<-rbind(pred2,x)
	}
pred2<-data.frame(pred2,ws)

pdf(file="pred_scarcity.pdf",4,10)
par(mfrow=c(2,1))
plot(ws,pred1$mean,type="l",ylim=c(0,0.4),ylab="Probability of Cooperative Interaction",xlab="Water Scarcity")
for(i in c(1:nrow(pred1))){
segments(x0=pred1[i,6], y0=pred1[i,4],x1=pred1[i,6],y1=pred1[i,5],col="lightgray")
}
lines(ws,pred1$mean)
plot(ws,pred2$mean,type="l",ylim=c(0,0.4),ylab="Probability of Conflictual Interaction",xlab="Water Scarcity")
for(i in c(1:nrow(pred2))){
segments(x0=pred2[i,6], y0=pred2[i,4],x1=pred2[i,6],y1=pred2[i,5],col="lightgray")
}
lines(ws,pred2$mean)
dev.off()

mx<-max(data$mdevsd,na.rm=TRUE)
mn<-min(data$mdevsd,na.rm=TRUE)


z.out<-coop2
x.max<-setx(z.out,mdevsd=mx,hist_coop12=hiscoop, grow=gr,secondhalf=sh,exchange=ex,syr=sy,energy_diff=en)
x.min<-setx(z.out,mdevsd=mn,hist_coop12=hiscoop, grow=gr,secondhalf=sh,exchange=ex,syr=sy,energy_diff=en)

s.out<-sim(z.out,x=x.min,x1=x.max,num=50000)
summary(s.out)

mean_coop<-mean(s.out$qi$fd)
quant_coop<-quantile(s.out$qi$fd,c(0.05,0.95))
low_coop<-quant_coop[1]
up_coop<-quant_coop[2]


z.out<-non2
x.max<-setx(z.out,mdevsd=mx,hist_non12=hisnon,grow=gr,secondhalf=sh,exchange=ex,syr=sy,energy_diff=en)
x.min<-setx(z.out,mdevsd=mn,hist_non12=hisnon,grow=gr,secondhalf=sh,exchange=ex,syr=sy,energy_diff=en)
s.out<-sim(z.out,x=x.min,x1=x.max,num=50000)
summary(s.out)

mean_non<-mean(s.out$qi$fd)
quant_non<-quantile(s.out$qi$fd,c(0.05,0.95))
low_non<-quant_non[1]
up_non<-quant_non[2]

mean<-cbind(mean_non,mean_coop)

#Plot of the first differences for the full basin sample

pdf(file="fdall.pdf",7,5)
plot(x = mean,y = c(1,2), xlim=c(-0.05,0.4),ylim=c(0.5,2.5),yaxt="n",xlab="predicted substantive effect of water scarcity (min to max) on probability of interaction", ylab="type of interaction")
axis(2,1:2,labels=c("conflictual","cooperative"))

abline(v=0)
segments(x0=low_coop,y0=2,x1=up_coop,y1=2,lty=2,col="black")
segments(x0=low_coop,y0=1.97,x1=low_coop,y1=2.03,col="black")
segments(x0=up_coop,y0=1.97,x1=up_coop,y1=2.03,col="black")

segments(x0=low_non,y0=1,x1=up_non,y1=1,col="gray")
segments(x0=low_non,y0=0.97,x1=low_non,y1=1.03,col="gray")
segments(x0=up_non,y0=0.97,x1=up_non,y1=1.03,col="gray")
dev.off()

##############################
## 		  					##
##		Split-basin  	    ##
##		  analyses			##
##							##
##############################


# Split Basin analyses
Syr<-subset(data, data$syr==1)
Amu<-subset(data, data$amu==1)


syr_coop<-clogit(Coop~mdevsd+hist_coop12+grow+secondhalf+ strata(dyad), data=Syr)
syr_non<-clogit(Non~mdevsd+hist_non12+grow+secondhalf+ strata(dyad), data=Syr)

amu_coop<-clogit(Coop~mdevsd+hist_coop12+grow+secondhalf+ strata(dyad), data=Amu)
amu_non<-clogit(Non~mdevsd+hist_non12+grow+secondhalf+ strata(dyad), data=Amu)


##############################
## 		  					##
##	  Robustness Checks  	##
##							##
##############################


### Using pre-2000 data

coop1<-zelig(Coop~mdevpresd, data=data, model="logit.gee",id="dyad")
coop2<-zelig(Coop~mdevpresd+hist_coop12+grow+secondhalf+ energy_diff+exchange + syr , data= data,model="logit.gee",id="dyad")
coop3<-clogit(Coop~mdevpresd+strata(dyad), data=data)
coop4<-clogit(Coop~mdevpresd+hist_coop12+grow+secondhalf+energy_diff + strata(dyad), data=data)

non1<-zelig(Non~mdevpresd, data=data, model="logit.gee",id="dyad")
non2<-zelig(Non~mdevpresd+hist_non12+grow+secondhalf+ energy_diff+ exchange + syr , data=data,model="logit.gee",id="dyad")
non3<-clogit(Non~mdevpresd+strata(dyad), data=data)
non4<-clogit(Non~mdevpresd+hist_non12+grow+secondhalf+ energy_diff+strata(dyad), data=data)


### Using indicator for extreme scarcity

coop1<-zelig(Coop~mdevex, data=data, model="logit.gee",id="dyad")
coop2<-zelig(Coop~mdevex+hist_coop12+grow+secondhalf + energy_diff+ exchange + syr, data=data,model="logit.gee",id="dyad")
coop3<-clogit(Coop~mdevex+strata(dyad), data=data)
coop4<-clogit(Coop~mdevex+hist_coop12+grow+secondhalf + energy_diff+strata(dyad), data=data)

non1<-zelig(Non~mdevex, data=data, model="logit.gee", id="dyad")
non2<-zelig(Non~mdevex+hist_non12+grow+secondhalf + energy_diff+ exchange + syr + energy_diff, data=data,model="logit.gee",id="dyad")
non3<-clogit(Non~mdevex+strata(dyad), data=data)
non4<-clogit(Non~mdevex+hist_non12+grow+secondhalf + energy_diff+strata(dyad), data=data)


### Including lagged scarcity

coop1<-zelig(Coop~mdevsd+mdevsdlag, data=data, model="logit.gee",id="dyad")
coop2<-zelig(Coop~mdevsd+mdevsdlag+hist_coop12+grow+secondhalf+ energy_diff+ exchange + syr, data=data,model="logit.gee",id="dyad")
coop3<-clogit(Coop~mdevsd+mdevsdlag+strata(dyad), data=data)
coop4<-clogit(Coop~mdevsd+mdevsdlag+hist_coop12+grow+secondhalf+ energy_diff+strata(dyad), data=data)

non1<-zelig(Non~mdevsd+mdevsdlag, data=data, model="logit",id="dyad")
non2<-zelig(Non~mdevsd+mdevsdlag+hist_non12+grow+secondhalf+ energy_diff+ exchange + syr , data=data,model="logit.gee",id="dyad")
non3<-clogit(Non~mdevsd+mdevsdlag+strata(dyad), data=data)
non4<-clogit(Non~mdevsd+mdevsdlag+hist_non12+grow+secondhalf+ energy_diff+strata(dyad), data=data)


### Using 3-month scarcity

coop1<-zelig(Coop~m3devsd, data=data, model="logit.gee",id="dyad")
coop2<-zelig(Coop~m3devsd+hist_coop12+grow+secondhalf+ energy_diff+exchange + syr , data=data,model="logit.gee",id="dyad")
coop3<-clogit(Coop~m3devsd+strata(dyad), data=data)
coop4<-clogit(Coop~m3devsd+hist_coop12+grow+secondhalf+energy_diff + strata(dyad), data=data)

non1<-zelig(Non~m3devsd, data=data, model="logit.gee",id="dyad")
non2<-zelig(Non~m3devsd+hist_non12+grow+secondhalf+ energy_diff+exchange + syr, data=data,model="logit.gee",id="dyad")
non3<-clogit(Non~m3devsd+strata(dyad), data=data)
non4<-clogit(Non~m3devsd+hist_non12+grow+secondhalf+energy_diff + strata(dyad), data=data)


### Including energy difference

coop1<-zelig(Coop~mdevsd, data=data, model="logit.gee",id="dyad")
coop2<-zelig(Coop~mdevsd+hist_coop12+grow+secondhalf+energy_diff+exchange + syr, data=data,model="logit.gee",id="dyad")
coop3<-clogit(Coop~mdevsd + strata(dyad), data=data)
coop4<-clogit(Coop~mdevsd+hist_coop12+grow+secondhalf+energy_diff+strata(dyad),data=data)

non1<-zelig(Non~mdevsd, data=data, model="logit.gee",id="dyad")
non2<-zelig(Non~mdevsd +hist_non12+grow+secondhalf + energy_diff+exchange + syr, data=data,model="logit.gee",id="dyad")
non3<-clogit(Non~mdevsd + strata(dyad), data=data)
non4<-clogit(Non~mdevsd+hist_non12+grow+secondhalf+energy_diff + strata(dyad),data=data)


### Including contiguity

coop1<-zelig(Coop~mdevsd+contig, data=data, model="logit.gee",id="dyad")
coop2<-zelig(Coop~mdevsd+contig+hist_coop12+ grow + secondhalf +  syr , data=data,model="logit.gee",id="dyad")

non1<-zelig(Non~mdevsd+contig, data=data, model="logit.gee",id="dyad")
non2<-zelig(Non~mdevsd+contig+hist_non12+ grow + secondhalf +  syr , data=data,model="logit.gee",id="dyad")


### Excluding failed negotiations from conflictual interactions

coop1<-zelig(Coop~mdevsd, data=nobreak, model="logit.gee",id="dyad")
coop2<-zelig(Coop~mdevsd+hist_coop12+grow+secondhalf + energy_diff+ exchange + syr, data=nobreak,model="logit.gee",id="dyad")
coop3<-clogit(Coop~mdevsd+strata(dyad), data=nobreak)
coop4<-clogit(Coop~mdevsd+hist_coop12+grow+secondhalf + energy_diff+strata(dyad), data=nobreak)

non1<-zelig(Non~mdevsd, data=nobreak, model="logit.gee", id="dyad")
non2<-zelig(Non~mdevsd+hist_non12+grow+secondhalf + energy_diff+ exchange + syr + energy_diff, data=nobreak,model="logit.gee",id="dyad")
non3<-clogit(Non~mdevsd+strata(dyad), data=nobreak)
non4<-clogit(Non~mdevsd+hist_non12+grow+secondhalf + energy_diff+strata(dyad), data=nobreak)


### Only water-related events

coop1<-zelig(Coop~mdevsd, data=water, model="logit.gee",id="dyad")
coop2<-zelig(Coop~mdevsd+hist_coop12 +grow+secondhalf+exchange + syr, data=water,model="logit.gee",id="dyad")
coop3<-clogit(Coop~mdevsd + strata(dyad), data=water)
coop4<-clogit(Coop~mdevsd+hist_coop12+grow+secondhalf+strata(dyad),data=water)

non1<-zelig(Non~mdevsd, data=water, model="logit.gee",id="dyad")
non2<-zelig(Non~mdevsd +hist_non12+grow+secondhalf +exchange + syr, data=water,model="logit.gee",id="dyad")
non3<-clogit(Non~mdevsd + strata(dyad), data=water)
non4<-clogit(Non~mdevsd+hist_non12+grow+secondhalf + strata(dyad),data=water)

### Multinomial logit models

data$drop<-rep(0, nrow(data))
for (i in c(2:nrow(data))){
	j<-i-1
	k<-ncol(data)
	data$drop[i]<-ifelse((data$dyad[i]==data$dyad[j] & data$Date[i] == data$Date[j]),1,0)
}

multi<-subset(data,data$drop==0)

#three categories
threecat<-zelig(as.factor(cat3)~mdevsd, model="mlogit",data=multi)
summary(threecat)
threecat2<-zelig(as.factor(cat3)~mdevsd+hist_coop12+hist_non12, model="mlogit",data=multi)
summary(threecat2)

#four categories
fourcatone<-zelig(as.factor(cat4)~mdevsd,model=
"mlogit",data=multi)
summary(fourcatone)
fourcattwo<-zelig(as.factor(cat4)~mdevsd+hist_coop12+hist_non12,model=
"mlogit",data=multi)
summary(fourcattwo)


### Survival analysis 

modcoop1<-coxph(Surv(start,stop,fail)~mdevsd ,data=survnon)
summary(modcoop1)
modnon1<-coxph(Surv(start,stop,fail)~mdevsd,data=survcoop)
summary(modnon1)

pdf(file="survival.pdf",8,8)
par(mfrow=c(2,1))
attach(survnon)
d<-data.frame(mdevsd=c(min(mdevsd,na.rm=TRUE),max(mdevsd,na.rm=TRUE)))
detach(survnon)
plot(survfit(modcoop1,newdata=d),lty=c(1,2),xlim=c(0,12),mark.time=FALSE,main = "survival function for conflictual state",xlab="months",ylab="proportion without cooperative interaction")
legend("topright",legend=c("max. water scarcity","min water scarcity"),lty=c(1,2))
attach(survcoop)
d<-data.frame(mdevsd=c(min(mdevsd,na.rm=TRUE),max(mdevsd,na.rm=TRUE)))
detach(survcoop)
plot(survfit(modnon1,newdata=d),lty=c(1,2),xlim=c(0,12),mark.time=FALSE, main="survival function for cooperative state",xlab="months",ylab ="proportion without conflictual inferaction")
legend("topright",legend=c("max. water scarcity","min. water scarcity"),lty=c(1,2))

dev.off()